1 JB changes

Removed DESeq part and added some file saving for CSVs

# Load libraries
library(DESeq2)
library(edgeR)
library(limma)
library(variancePartition)
library(foreach)
# library(clusterProfiler)
library(org.Hs.eg.db)
# library(ggh4x)
library(magrittr)
library(enrichplot)
library(patchwork)
library(fgsea)
library(doParallel)
library(tidyverse)

# Define main path
main_path <- file.path("/home/jbrenton/ASAP_Hardy_post_align_files/Final_Dataset/Hardy_Limma_Analysis_Final_version_only/")

# main_path <-  here::here()

# Load helper functions
source(file.path(main_path, "R/hf_graph_and_themes.R"))
source(file.path(main_path, "R/hf_reports.R"))

# Set options
options(dplyr.summarise.inform = FALSE)
options(lifecycle_verbosity = "warning")
knitr::opts_chunk$set(echo = F, warning = T, message = T, 
                      out.width="100%", fig.align = "center", dpi = 150,
                      fig.width = 7.2, fig.height = 7.2)

## Custom sciRmdTheme with increased width.
## Please ask guillermorocamora@gmail.com for more details.
## Download from devtools::install_github("https://github.com/guillermo1996/grpSciRmdTheme")
if(require(grpSciRmdTheme)){
  grpSciRmdTheme::set.theme(
    theme = "default",
    color = NULL,
    header.sticky = FALSE,
    list.group.icon = "arrow",
    font.family = "Arial",
    font.color = "black",
    header.color = "darkblue",
    figure.perc = 100
  )
}

## Custom styles for the output html
cat('
<style type="text/css">
.dataTables_scrollHeadInner{
  width:100% !important;
}
.dataTables_scrollHeadInner table{
  width:100% !important;
}
/*
.code-folding-btn {
  display: none;
}
*/
h3, .h3 {
  font-size: 22px!important;
}
h4, .h4 {
  font-size: 18px!important;
}
h5, .h5 {
  font-size: 16px!important;
}
body{
  font-size: 13px;
}

.tocify-subheader {
    text-indent: 15px;
    display: none;
    font-size: 12px;
}

.break5 .nav-tabs, .break4 .nav-tabs, .break3 .nav-tabs, .break2 .nav-tabs{
  display: grid;
  text-align: center;
  font-size: 14px;
  font-weight: bold;
}
.break5 .nav-tabs:before, .break4 .nav-tabs:before, .break3 .nav-tabs:before, .break2 .nav-tabs:before,
.break5 .nav-tabs:after, .break4 .nav-tabs:after, .break3 .nav-tabs:after, .break2 .nav-tabs:after {
  content: unset;
}

.break5 .nav-tabs{
  grid-template-columns: repeat(5, 1fr);
}
.break4 .nav-tabs{
  grid-template-columns: repeat(4, 1fr);
}
.break3 .nav-tabs{
  grid-template-columns: repeat(3, 1fr);
}
.break2 .nav-tabs{
  grid-template-columns: repeat(2, 1fr);
}

.breakB5 .nav-tabs{
  grid-template-columns: repeat(5, 1fr);
  font-size: 12px;
}
.breakB4 .nav-tabs{
  grid-template-columns: repeat(4, 1fr);
  font-size: 12px;
}
.breakB3 .nav-tabs{
  grid-template-columns: repeat(3, 1fr);
  font-size: 12px;
}
.breakB2 .nav-tabs{
  grid-template-columns: repeat(2, 1fr);
  font-size: 12px;
}

.stop-break .nav-tabs {
  display: flex;
  font-size: 13px;
  font-weight: normal;
}
.stop-break .nav-tabs:before,
.stop-break .nav-tabs:after {
  content: unset;
}
</style>')

2 Background

This report includes a differential gene expression analysis on the Hardy bulk dataset with limma package, using voom and duplicateCorrelation. It is created as a template and not as a report with the final results.

2.1 Load previous data

The main path to run this report on the RytenLab server is at:

/home/grocamora/RytenLab-Research/18-DGE_limma_dream/

This report is meant to be run after after:

Covariate identification script: located at Scripts/01-Hardy_covariate_identification.R
Covariate identification report: located at RMarkdown/02-Hardy_covariate_identification.R
PCA Analysis for outliers: located at RMarkdown/03-PC_Analysis_Outliers.Rmd

The most relevant data to load prior to the analysis are the metadata, the salmon quantifications, a list of blacklisted genes and the outliers to remove:

# Input paths
main_path <- "/home/jbrenton/ASAP_Hardy_post_align_files/Final_Dataset/Hardy_Limma_Analysis_Final_version_only/"
results_path <- file.path(main_path, "results")

txi.salmon_path <- file.path(main_path, "data/Hardy/txi.salmon.rds")
sample_metadata_path <- file.path(main_path, "metadata/Hardy/updated_metadata_w_PMI.csv")
metadata_path <- file.path(results_path, "all_brain_areas_covariates.rds")
sample_outliers_path <- file.path(results_path, "sample_outliers16.rds")

blacklist_genes <- file.path(main_path, "data/blist_genes_no_vers.csv")
gene_map_path <- file.path(main_path, "data/gencode_txid_to_geneid.txt")

# Output paths
hardy_dge_results_path <- file.path(results_path, "dge_variables/")

deg_counts_output_path <- file.path(hardy_dge_results_path, "PD_Collapsed_dge_counts.rds")
dge_all_path <- file.path(hardy_dge_results_path, "PD_Collapsed_dge_all.rds")

dir.create(hardy_dge_results_path, showWarnings = F, recursive = F)

# Load data
txi.salmon <- readRDS(txi.salmon_path)
cts <- txi.salmon$counts

metadata <- loadMetadataHardy(metadata_path, sample_metadata_path)
sample_outliers <- readRDS(sample_outliers_path)
bxp_outliers <- sample_outliers$bxp_id_full[sample_outliers$outlier]

# Load blacklisted genes
blacklist_genes <- vroom::vroom(blacklist_genes, col_names = T, delim = ",")
blacklist_genes <- as.vector(unlist(blacklist_genes))

gencode_txid_to_geneid <- vroom::vroom(gene_map_path) %>%
  `colnames<-`(c("tx_id", "gene_id","gene_name", "description")) %>%
  dplyr::mutate(tx_id = sub("\\..+", "", tx_id),
                gene_id = sub("\\..+", "", gene_id))

# Recode PD to combine PD groups
metadata %<>% mutate(pd=case_when(pd=="PD" ~ "PD",
                                 pd=="PDD" ~ "PD",
                                 pd=="control" ~ "control"))

The covariates employed to model the gene expression data are extracted from the Covariate identification report previously mentioned and are the following:

GC_NC_40_59, PCT_UTR_BASES, gender, Total_Sequences, INTRONIC_BASES and tss_up_1kb_tag_pct.

The experimental condition of interest in this report will be the type of PD of the samples: controls, PD without dementia (PD) and PD with dementia (PDD). We will also study two different scenarios for the brain regions: i) a collapsed brain region approach and ii) a per brain region analysis.

experimental_condition = "pd"
experimental_region = "brain_region"

It is also required to fill in certain parameters that will affect the pipeline and the results:

gene_filter_method: different possible gene expression thresholds can be tested by modifying this variable. These are explained in the preprocessing section. Only possible values are DEGreport, edgeR and cpm. If not a valid value is provided, no gene threshold is applied.

# Gene filtering approach
gene_filter_method = "DEGreport"

if(!gene_filter_method %in% c("DEGreport", "edgeR", "cpm")){
  gene_filter_method = "None"
  warning("not a valid gene filter method was supplied. Defaults to no threshold applied.")
}

Changelog

v1.6

Updated to R v4.3.0

v1.5

Fixed covariate identification mistake.

v1.4

Updated to final decision on outliers and gene expression threshold.
Removed Gene Enrichment Analyses from this report.

v1.3

Added GO, KEGG and REACTOME gene enrichment analyses.
Updated default report to include 16 outliers with the DGEreport::degFilter threshold.

v1.2

Merged the two limma pipelines to combine diagnosis and brain region. Common practice in the literature.
Added DEGreport filter to the gene expression threhsolds.
Added a summary table of all results.

v1.1

Removed the dream pipeline.
Added the limma/voom pipeline for all tissues combined and by brain region.
Added DESeq2 analysis for reference.
Added additional study about the gene expression threshold.

v1.0

Initial report

3 DGE with `limma/voom` and `duplicateCorrelation`

A complete tutorial on how to use limma can be found in their documentation page. On its core, limma fits a linear model to assess differential gene expression across different comparisons (or contrasts).

Given that we are comparing across different tissues with samples from the same individuals, we need to treat the sample donor as a random effect. In limma, this is handled by the duplicateCorrelation function, which estimates a correlation between replicates from the same individual and gene expressions. In this approach, a single correlation value is employed genome-wide, assuming that the correlation is the same across all genes. Other alternatives, such as the dream package, estimate a per-gene correlation, but their study is beyond the scope of this report.

Following the instructions in the official documentation and in this tutorial about design matrices, the two factors of interests (PD status and brain region) are grouped into a single factor. Thus, only one model is fitted and all comparisons are extracted from it.

3.1 Preprocessing

The first step is to generate a dataframe with the covariates of interests (limma_covs) and a dataframe with the covariates plus the grouping factors (limma_covs_all) which includes the experimental condition, the brain region and the individual ID. Outliers are also removed following the criteria explained in PCA Analysis for outliers report: outliers detected using PCA Z-scores and WGCNA sample connectivity Z-scores. It is important to ensure that the order of samples in the counts object is the same as the order of samples in the metadata table.

# Ensure sample order
metadata <- metadata[sapply(metadata$bxp_id_full, function(x) grep(x, colnames(cts))), ]
stopifnot(all(colnames(cts) == metadata$bxp_id_full))

# Covariate + grouping dataframe
limma_covs_all <- metadata %>%
  dplyr::select(sample_id, bxp_id_full, individual_id = case_id,
                condition = all_of(experimental_condition),
                region = all_of(experimental_region),
                all_of(covariates)) %>%
  # Filter the outliers
  dplyr::filter(!bxp_id_full %in% bxp_outliers) %>%
  # Generate the grouping factor as "condition_region"
  dplyr::mutate(group = as.factor(paste0(condition, "_", region))) %>%
  dplyr::relocate(individual_id, group) %>%
  # Convert all characters to factor
  dplyr::mutate(across(where(is.character), ~ as.factor(.))) %>%
  # Scale the covariates again after the removal of outliers
  dplyr::mutate(across(where(is.numeric), ~ scale(., center = T, scale = T))) %>%
  tibble::column_to_rownames("bxp_id_full")

# Covariates dataframe
limma_covs <- limma_covs_all %>% dplyr::select(all_of(covariates))

Next, we generate the DGE object and apply the following preprocessing steps:

Remove outliers.
Remove blacklisted genes.
Remove low expressed genes: there are currently four different options:
- DEGreport: employs the function DEGreport::degFilter, which requires at least one count in all samples of at least one group.
- edgeR: employs the function edgeR::filterByExpr which, in this situation, requires a minimum cpm of 0.11 in at least 19 samples.
- cpm: requires a minimum cpm of one in at least 50% of the samples.
- None: no gene expression filter is applied.
Calculate the normalization factors.

# The author recommends to use "keep.lib.sizes = FALSE" after gene filtering,
# but before "calcNormFactors". Source:
# https://support.bioconductor.org/p/9136787/

# Remove outliers
cts_clean <- cts[, !colnames(cts) %in% bxp_outliers]

# Generate the DGE object
dge0 <- edgeR::DGEList(cts_clean, group = limma_covs_all$group)

# Remove the blacklisted genes
dge1 <- dge0[!rownames(dge0) %in% blacklist_genes, , keep.lib.sizes = F]

# Remove low-expressed genes. We are currently testing different approaches to
# filtering the low-expressed genes. Three valid approaches are currently
# employed:
#
# i) edgeR: use the function "edgeR::filterByExpr". It currently requires a cpm
# > 0.11 in at least 19 samples.
#
# ii) cpm: the same filter applied to covariate identification and PC analysis.
# It requires a cpm or at least 1 in 50% of the samples.
#
# iii) DEGreport: a previously used filter by the group. It requires at least
# one count across all samples in one of the groups.

if(gene_filter_method == "DEGreport"){
  filter_count <- DEGreport::degFilter(counts = dge1$counts,
                                       metadata = limma_covs_all,
                                       group = "group",
                                       min = 1, # All samples in group must have more than expr > 0
                                       minreads = 0)
  isexpr <- rownames(dge1$counts) %in% rownames(filter_count)
}else if(gene_filter_method == "edgeR"){
  isexpr <- edgeR::filterByExpr(dge1, group = limma_covs_all$group)
}else if(gene_filter_method == "cpm"){
  isexpr <- rowSums(edgeR::cpm(dge1) > 1) >= 0.5*ncol(dge1)
}else{
  isexpr <- rep(T, nrow(dge1$counts))
}

dge2 <- dge1[isexpr, , keep.lib.sizes = F]

# Calculate normalization factors - adds that column to the samples list
dge3 <- edgeR::calcNormFactors(dge2)

# Final DGE object
dge <- dge3

# Clear variables
rm(cts_clean, dge0, dge1, dge2, dge3)

message(paste0("Employed gene expression filter: ", gene_filter_method))

{bg.success} ## Employed gene expression filter: DEGreport

Finally, we generate the design object required for the limma pipeline. We also specify the formula that the linear model fit will use, and it must be shaped as ~ 0 + experimental_variable + covariates. The use of intercept is advised to allow for all different comparisons. The formula employed in the analysis is the following:

Formula: ~ 0 + group + gender + GC_NC_40_59 + tss_up_1kb_tag_pct + PCT_UTR_BASES + Total_Sequences + INTRONIC_BASES

limma_formula <- reformulate(c(0, "group", limma_covs %>% dplyr::select(where(is.factor), where(is.numeric)) %>% colnames()))
limma_factor_formula <- reformulate(c(0, "group", limma_covs %>% dplyr::select(where(is.factor)) %>% colnames()))
limma_num_formula <- reformulate(limma_covs %>% dplyr::select(where(is.numeric)) %>% colnames())

limma_mod <- model.matrix(limma_factor_formula, data = limma_covs_all) %>%
  `colnames<-`(stringr::str_remove_all(colnames(.), "^group"))
limma_design <- data.frame(limma_mod, limma_covs %>% dplyr::select(where(is.numeric)))

3.2 Voom with `duplicateCorrelation`

We use voom to obtain the weights of each gene based on the residuals of a linear fit (more information). Since we also want to consider the correlation between samples from the same individual, we employ duplicateCorrelation. Both functions are applied twice based on recommendations from the developers.

By default, the following code block stores and load from memory the vobj and dupcor variables. Please remove files from disk or modify the results folder to run the analysis again.

vobj_path <- file.path(hardy_dge_results_path, "PD_Collapsed_vobj_limma.rds")
dupcor_path <- file.path(hardy_dge_results_path, "PD_Collapsed_dupcor_limma.rds")

if(!file_test("-f", vobj_path, dupcor_path)){
  vobj_tmp <- limma::voom(dge, limma_design, save.plot = T)
  dupcor_tmp <- limma::duplicateCorrelation(vobj_tmp, limma_design, block = limma_covs_all$individual_id)

  vobj <- limma::voom(dge, limma_design, save.plot = T,
                      block = limma_covs_all$individual_id, correlation = dupcor_tmp$consensus)
  dupcor <- limma::duplicateCorrelation(vobj, limma_design, block = limma_covs_all$individual_id)


  vobj %>% saveRDS(vobj_path)
  dupcor %>% saveRDS(dupcor_path)
}else{
  vobj <- readRDS(vobj_path)
  dupcor <- readRDS(dupcor_path)
}

plotVobj(vobj)

3.3 Fit model

The next step is to fit the model for each gene. We employ the function lmFit followed by eBayes to estimate the log-odds of differential expression and t-statistic.

fit_path <- file.path(hardy_dge_results_path, "PD_Collapsed_fit_limma.rds")
if(!file_test("-f", fit_path)){
  fit <- limma::lmFit(vobj, limma_design, block = limma_covs_all$individual_id, correlation = dupcor$consensus)
  fit <- limma::eBayes(fit)

  fit %>% saveRDS(fit_path)
}else{
  fit <- readRDS(fit_path)
}

3.4 Make constrasts

One of the most convenient aspects of the limma pipeline is the ability to accommodate arbitrary experimental complexity. The function limma::contrasts.fit allows for comparisons between groups to be made after the fit. To generate these contrasts, we need to define a contrast matrix specifing the groups to compare. It is important to notice that each positive contribution must add to 1 and each negative contribution must add to -1 (otherwise, the log fold-changes are scaled). The reference group must have a negative contribution.

contrasts_fit <- limma::makeContrasts(
  "Collapse-PD-Control" = (PD_ACG + PD_IPL + PD_MFG + PD_MTG -
                           control_ACG - control_IPL - control_MFG - control_MTG)/4,
    "ACG-PD-Control" = PD_ACG - control_ACG,
  "IPL-PD-Control" = PD_IPL - control_IPL,
   "MFG-PD-Control" = PD_MFG - control_MFG,
  "MTG-PD-Control" = PD_MTG - control_MTG,
  levels = limma_design)

# Fit by the contrasts
fit_con <- limma::eBayes(limma::contrasts.fit(fit, contrasts_fit))

## Warning in limma::contrasts.fit(fit, contrasts_fit): row names of contrasts
## don't match col names of coefficients

In our scenario, a total of 15 different contrasts were made: 3 different case vs. control studies (i.e. PD vs control, PDD vs control and PDD vs PD) in 5 different tissue combinations (collapsed and by brain region). The following table summarizes the contrasts:

4 Results: Collapsed brain region

In this section, we will compare the results of combining all brain regions. The generated contrast table is drawn below, where each coefficient is set to \(1/\text{number of brain regions}\). As mentioned previously, the sign of the coefficient is important, leaving the negative contributions to the reference group.

To extract the results from the fit object, we employ a custom function that internally calls limma::topTable for a given coefficient (i.e. the name for the contrasts of interest in the contrasts matrix):

dge_collapsed_i <- extractResultsDGE(fit_con, coef = "Collapse-PD-Control")

dge_collapsed <- dplyr::bind_rows(
  dge_collapsed_i %>% dplyr::mutate(tissue = "Collapse", test = "PD-Control"))

The number of differentially expressed genes obtained for each comparison is summarized in the following table:

##   Down NotSig     Up 
##    123  22786     62

4.1 Volcano plots

Volcano plots allow for a more global view of the DGE results by showing statistical significance vs magnitude of change - i.e. log transformed p-values in the Y-axis and the log2 fold-change in the X-axis. Our requirement for significance is \(p-value <= 0.05\), drawn as a red line. Every gene above the mentioned line is considered differentially expressed. The sign of the log2 fold-change specifies whether it is upregulated or downregulated. Additionally, the top most significant differentially expressed genes are named in each test.

PD vs Control

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_text_repel()`).

4.2 Top significant genes

Here we represent the 10 most significant differentially expressed genes along with their magnitude of change and their statistical significance - i.e. log transformed p-values as the fill colour and log2 fold-change in the Y-axis.

PD vs Control

4.3 Mean-difference plots

The mean-difference plots are another useful tool to visualize the DGE results. Instead of showing the statistical significance (as the Volcano plot), they represent the magnitude of change against the average log-CPM expression - i.e. average log-CPM in the x-Axis and log2 fold-change in the X axis. Significant differentially expressed genes are represented both by their color and their shape:

PD vs Control

5 Results: By brain region

# Extract the results from the main fit.
tissue_contrasts <- colnames(contrasts_fit)[-grep("Collapse", colnames(contrasts_fit))]
dge_tissues <- lapply(tissue_contrasts, function(contrast_name){
  tissue <- stringr::str_split_fixed(contrast_name, "-", n = 2)[1]
  test <- stringr::str_split_fixed(contrast_name, "-", n = 2)[2]

  extractResultsDGE(fit_con, coef = contrast_name) %>% dplyr::mutate(tissue = tissue, test = test)
}) %>% dplyr::bind_rows()

In this section, we will compare the results for the different brain regions.

Expand to see contrast table.

Expand to see number of differentially expressed genes.

## Warning: 'xfun::attr()' is deprecated.
## Use 'xfun::attr2()' instead.
## See help("Deprecated")

Tissue - Test	Downregulated	Unchanged	Upregulated
ACG
PD-Control	48	22909	14
IPL
PD-Control	17	22952	2
MFG
PD-Control	33	22932	6
MTG
PD-Control	49	22870	52

5.1 Volcano plots

ACG

IPL

MFG

MTG

5.2 Top Significant genes

ACG

IPL

MFG

MTG

5.3 Mean-difference plots

ACG

IPL

MFG

MTG

6 Summary of results

In this section, a brief summary of the results will be presented.

6.1 Number of differentially expressed genes

First, we can study the number of differentially expressed genes, regardless of their magnitude or direction of change:

Expand to see proportion of upregulated/downregulated

Each tile contains the proportion of upregulated/downregulated genes. Only proportions higher than 30% are written in the tiles for clarity:

PD vs Control

## Warning: 'xfun::attr()' is deprecated.
## Use 'xfun::attr2()' instead.
## See help("Deprecated")

Tissue	Regulated	Number of DEG - Limma
Collapse
Collapse	Significantly Downregulated	123
Collapse	Significantly Upregulated	62
ACG
ACG	Significantly Downregulated	48
ACG	Significantly Upregulated	14
IPL
IPL	Significantly Downregulated	17
IPL	Significantly Upregulated	2
MFG
MFG	Significantly Downregulated	33
MFG	Significantly Upregulated	6
MTG
MTG	Significantly Downregulated	49
MTG	Significantly Upregulated	52

6.2 Common DEGs between brain regions

Expand to see results of including the collapsed approach.

7 Output results to csv

8 Session Info

Click to Show

## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value
##  version  R version 4.3.0 (2023-04-21)
##  os       Ubuntu 20.04.6 LTS
##  system   x86_64, linux-gnu
##  ui       X11
##  language (EN)
##  collate  en_GB.UTF-8
##  ctype    en_GB.UTF-8
##  tz       Etc/UTC
##  date     2025-09-30
##  pandoc   3.1.1 @ /usr/lib/rstudio-server/bin/quarto/bin/tools/ (via rmarkdown)
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  package              * version    date (UTC) lib source
##  abind                  1.4-8      2024-09-12 [1] CRAN (R 4.3.0)
##  AnnotationDbi        * 1.64.1     2023-11-03 [1] Bioconductor
##  aod                    1.3.3      2023-12-13 [1] CRAN (R 4.3.0)
##  ape                    5.8-1      2024-12-16 [1] CRAN (R 4.3.0)
##  aplot                  0.2.4      2024-12-17 [1] CRAN (R 4.3.0)
##  backports              1.5.0      2024-05-23 [1] CRAN (R 4.3.0)
##  Biobase              * 2.62.0     2023-10-24 [1] Bioconductor
##  BiocGenerics         * 0.48.1     2023-11-01 [1] Bioconductor
##  BiocParallel         * 1.36.0     2023-10-24 [1] Bioconductor
##  Biostrings             2.70.3     2024-03-13 [1] Bioconductor 3.18 (R 4.3.0)
##  bit                    4.6.0      2025-03-06 [1] CRAN (R 4.3.0)
##  bit64                  4.6.0-1    2025-01-16 [1] CRAN (R 4.3.0)
##  bitops                 1.0-9      2024-10-03 [1] CRAN (R 4.3.0)
##  blob                   1.2.4      2023-03-17 [1] CRAN (R 4.3.0)
##  bookdown               0.42       2025-01-07 [1] CRAN (R 4.3.0)
##  boot                   1.3-31     2024-08-28 [1] CRAN (R 4.3.0)
##  broom                  1.0.7      2024-09-26 [1] CRAN (R 4.3.0)
##  bslib                  0.8.0      2024-07-29 [1] CRAN (R 4.3.0)
##  cachem                 1.1.0      2024-05-16 [1] CRAN (R 4.3.0)
##  caTools                1.18.3     2024-09-04 [1] CRAN (R 4.3.0)
##  circlize               0.4.16     2024-02-20 [1] CRAN (R 4.3.0)
##  cli                    3.6.5      2025-04-23 [1] CRAN (R 4.3.0)
##  clue                   0.3-66     2024-11-13 [1] CRAN (R 4.3.0)
##  cluster                2.1.8.1    2025-03-12 [1] CRAN (R 4.3.0)
##  codetools              0.2-20     2024-03-31 [1] CRAN (R 4.3.0)
##  colorspace             2.1-1      2024-07-26 [1] CRAN (R 4.3.0)
##  ComplexHeatmap         2.18.0     2023-10-24 [1] Bioconductor
##  ConsensusClusterPlus   1.66.0     2023-10-24 [1] Bioconductor
##  corpcor                1.6.10     2021-09-16 [1] CRAN (R 4.3.0)
##  cowplot                1.1.3      2024-01-22 [1] CRAN (R 4.3.0)
##  crayon                 1.5.3      2024-06-20 [1] CRAN (R 4.3.0)
##  crosstalk              1.2.1      2023-11-23 [1] CRAN (R 4.3.0)
##  data.table             1.17.6     2025-06-17 [1] CRAN (R 4.3.0)
##  DBI                    1.2.3      2024-06-02 [1] CRAN (R 4.3.0)
##  DEGreport              1.38.5     2023-12-06 [1] Bioconductor 3.18 (R 4.3.0)
##  DelayedArray           0.28.0     2023-10-24 [1] Bioconductor
##  DESeq2               * 1.42.1     2024-03-06 [1] Bioconductor 3.18 (R 4.3.0)
##  dichromat              2.0-0.1    2022-05-02 [1] CRAN (R 4.3.0)
##  digest                 0.6.37     2024-08-19 [1] CRAN (R 4.3.0)
##  doParallel           * 1.0.17     2022-02-07 [1] CRAN (R 4.3.0)
##  DOSE                   3.28.2     2023-12-10 [1] Bioconductor 3.18 (R 4.3.0)
##  dplyr                * 1.1.4      2023-11-17 [1] CRAN (R 4.3.0)
##  DT                     0.33       2024-04-04 [1] CRAN (R 4.3.0)
##  edgeR                * 4.0.16     2024-02-18 [1] Bioconductor 3.18 (R 4.3.0)
##  enrichplot           * 1.22.0     2023-10-24 [1] Bioconductor
##  EnvStats               3.0.0      2024-08-24 [1] CRAN (R 4.3.0)
##  evaluate               1.0.4      2025-06-18 [1] CRAN (R 4.3.0)
##  fANCOVA                0.6-1      2020-11-13 [1] CRAN (R 4.3.0)
##  farver                 2.1.2      2024-05-13 [1] CRAN (R 4.3.0)
##  fastmap                1.2.0      2024-05-15 [1] CRAN (R 4.3.0)
##  fastmatch              1.1-6      2024-12-23 [1] CRAN (R 4.3.0)
##  fgsea                * 1.28.0     2023-10-24 [1] Bioconductor
##  forcats              * 1.0.0      2023-01-29 [1] CRAN (R 4.3.0)
##  foreach              * 1.5.2      2022-02-02 [1] CRAN (R 4.3.0)
##  fs                     1.6.6      2025-04-12 [1] CRAN (R 4.3.0)
##  generics               0.1.4      2025-05-09 [1] CRAN (R 4.3.0)
##  GenomeInfoDb         * 1.38.8     2024-03-15 [1] Bioconductor 3.18 (R 4.3.0)
##  GenomeInfoDbData       1.2.11     2023-11-20 [1] Bioconductor
##  GenomicRanges        * 1.54.1     2023-10-29 [1] Bioconductor
##  GetoptLong             1.0.5      2020-12-15 [1] CRAN (R 4.3.0)
##  ggdendro               0.2.0      2024-02-23 [1] CRAN (R 4.3.0)
##  ggforce                0.4.2      2024-02-19 [1] CRAN (R 4.3.0)
##  ggfun                  0.1.8      2024-12-03 [1] CRAN (R 4.3.0)
##  ggplot2              * 4.0.0      2025-09-11 [1] CRAN (R 4.3.0)
##  ggplotify              0.1.2      2023-08-09 [1] CRAN (R 4.3.0)
##  ggraph                 2.2.1      2024-03-07 [1] CRAN (R 4.3.0)
##  ggrepel                0.9.6      2024-09-07 [1] CRAN (R 4.3.0)
##  ggtree                 3.10.1     2024-02-25 [1] Bioconductor 3.18 (R 4.3.0)
##  GlobalOptions          0.1.2      2020-06-10 [1] CRAN (R 4.3.0)
##  glue                   1.8.0      2024-09-30 [1] CRAN (R 4.3.0)
##  GO.db                  3.18.0     2023-11-23 [1] Bioconductor
##  GOSemSim               2.28.1     2024-01-17 [1] Bioconductor 3.18 (R 4.3.0)
##  gplots                 3.2.0      2024-10-05 [1] CRAN (R 4.3.0)
##  graphlayouts           1.2.1      2024-11-18 [1] CRAN (R 4.3.0)
##  gridExtra              2.3        2017-09-09 [1] CRAN (R 4.3.0)
##  gridGraphics           0.5-1      2020-12-13 [1] CRAN (R 4.3.0)
##  grpSciRmdTheme       * 0.4.0.1    2023-12-03 [1] Github (guillermo1996/grpSciRmdTheme@5f7645e)
##  gtable                 0.3.6      2024-10-25 [1] CRAN (R 4.3.0)
##  gtools                 3.9.5      2023-11-20 [1] CRAN (R 4.3.0)
##  HDO.db                 0.99.1     2023-11-23 [1] Bioconductor
##  hms                    1.1.3      2023-03-21 [1] CRAN (R 4.3.0)
##  htmltools              0.5.8.1    2024-04-04 [1] CRAN (R 4.3.0)
##  htmlwidgets            1.6.4      2023-12-06 [1] CRAN (R 4.3.0)
##  httr                   1.4.7      2023-08-15 [1] CRAN (R 4.3.0)
##  igraph                 2.1.4      2025-01-23 [1] CRAN (R 4.3.0)
##  IRanges              * 2.36.0     2023-10-24 [1] Bioconductor
##  iterators            * 1.0.14     2022-02-05 [1] CRAN (R 4.3.0)
##  jquerylib              0.1.4      2021-04-26 [1] CRAN (R 4.3.0)
##  jsonlite               2.0.0      2025-03-27 [1] CRAN (R 4.3.0)
##  kableExtra             1.4.0      2024-01-24 [1] CRAN (R 4.3.0)
##  KEGGREST               1.42.0     2023-10-24 [1] Bioconductor
##  KernSmooth             2.23-26    2025-01-01 [1] CRAN (R 4.3.0)
##  knitr                  1.49       2024-11-08 [1] CRAN (R 4.3.0)
##  labeling               0.4.3      2023-08-29 [1] CRAN (R 4.3.0)
##  lattice                0.22-7     2025-04-02 [1] CRAN (R 4.3.0)
##  lazyeval               0.2.2      2019-03-15 [1] CRAN (R 4.3.0)
##  lifecycle              1.0.4      2023-11-07 [1] CRAN (R 4.3.0)
##  limma                * 3.58.1     2023-10-31 [1] Bioconductor
##  lme4                   1.1-36     2025-01-11 [1] CRAN (R 4.3.0)
##  lmerTest               3.1-3      2020-10-23 [1] CRAN (R 4.3.0)
##  locfit                 1.5-9.10   2024-06-24 [1] CRAN (R 4.3.0)
##  logging                0.10-108   2019-07-14 [1] CRAN (R 4.3.0)
##  lubridate            * 1.9.4      2024-12-08 [1] CRAN (R 4.3.0)
##  magrittr             * 2.0.4      2025-09-12 [1] CRAN (R 4.3.0)
##  MASS                   7.3-60.0.1 2024-01-13 [1] CRAN (R 4.3.0)
##  Matrix                 1.6-5      2024-01-11 [1] CRAN (R 4.3.0)
##  MatrixGenerics       * 1.14.0     2023-10-24 [1] Bioconductor
##  matrixStats          * 1.5.0      2025-01-07 [1] CRAN (R 4.3.0)
##  memoise                2.0.1      2021-11-26 [1] CRAN (R 4.3.0)
##  minqa                  1.2.8      2024-08-17 [1] CRAN (R 4.3.0)
##  mnormt                 2.1.1      2022-09-26 [1] CRAN (R 4.3.0)
##  mvtnorm                1.3-3      2025-01-10 [1] CRAN (R 4.3.0)
##  nlme                   3.1-166    2024-08-14 [1] CRAN (R 4.3.0)
##  nloptr                 2.2.1      2025-03-17 [1] CRAN (R 4.3.0)
##  numDeriv               2016.8-1.1 2019-06-06 [1] CRAN (R 4.3.0)
##  org.Hs.eg.db         * 3.18.0     2023-11-23 [1] Bioconductor
##  patchwork            * 1.3.0      2024-09-16 [1] CRAN (R 4.3.0)
##  pbkrtest               0.5.3      2024-06-26 [1] CRAN (R 4.3.0)
##  pillar                 1.11.0     2025-07-04 [1] CRAN (R 4.3.0)
##  pkgconfig              2.0.3      2019-09-22 [1] CRAN (R 4.3.0)
##  plyr                   1.8.9      2023-10-02 [1] CRAN (R 4.3.0)
##  png                    0.1-8      2022-11-29 [1] CRAN (R 4.3.0)
##  polyclip               1.10-7     2024-07-23 [1] CRAN (R 4.3.0)
##  psych                  2.4.12     2024-12-23 [1] CRAN (R 4.3.0)
##  purrr                * 1.1.0      2025-07-10 [1] CRAN (R 4.3.0)
##  qvalue                 2.34.0     2023-10-24 [1] Bioconductor
##  R6                     2.6.1      2025-02-15 [1] CRAN (R 4.3.0)
##  rbibutils              2.3        2024-10-04 [1] CRAN (R 4.3.0)
##  RColorBrewer           1.1-3      2022-04-03 [1] CRAN (R 4.3.0)
##  Rcpp                   1.1.0      2025-07-02 [1] CRAN (R 4.3.0)
##  RCurl                  1.98-1.17  2025-03-22 [1] CRAN (R 4.3.0)
##  Rdpack                 2.6.4      2025-04-09 [1] CRAN (R 4.3.0)
##  readr                * 2.1.5      2024-01-10 [1] CRAN (R 4.3.0)
##  reformulas             0.4.0      2024-11-03 [1] CRAN (R 4.3.0)
##  remaCor                0.0.18     2024-02-08 [1] CRAN (R 4.3.0)
##  reshape                0.8.10     2025-06-19 [1] CRAN (R 4.3.0)
##  reshape2               1.4.4      2020-04-09 [1] CRAN (R 4.3.0)
##  RhpcBLASctl            0.23-42    2023-02-11 [1] CRAN (R 4.3.0)
##  rjson                  0.2.23     2024-09-16 [1] CRAN (R 4.3.0)
##  rlang                  1.1.6      2025-04-11 [1] CRAN (R 4.3.0)
##  rmarkdown              2.29       2024-11-04 [1] CRAN (R 4.3.0)
##  RSQLite                2.3.9      2024-12-03 [1] CRAN (R 4.3.0)
##  rstudioapi             0.17.1     2024-10-22 [1] CRAN (R 4.3.0)
##  S4Arrays               1.2.1      2024-03-04 [1] Bioconductor 3.18 (R 4.3.0)
##  S4Vectors            * 0.40.2     2023-11-23 [1] Bioconductor 3.18 (R 4.3.0)
##  S7                     0.2.0      2024-11-07 [1] CRAN (R 4.3.0)
##  sass                   0.4.9      2024-03-15 [1] CRAN (R 4.3.0)
##  scales               * 1.4.0      2025-04-24 [1] CRAN (R 4.3.0)
##  scatterpie             0.2.4      2024-08-28 [1] CRAN (R 4.3.0)
##  sessioninfo            1.2.2      2021-12-06 [1] CRAN (R 4.3.0)
##  shadowtext             0.1.4      2024-07-18 [1] CRAN (R 4.3.0)
##  shape                  1.4.6.1    2024-02-23 [1] CRAN (R 4.3.0)
##  SparseArray            1.2.4      2024-02-11 [1] Bioconductor 3.18 (R 4.3.0)
##  statmod                1.5.0      2023-01-06 [1] CRAN (R 4.3.0)
##  stringi                1.8.7      2025-03-27 [1] CRAN (R 4.3.0)
##  stringr              * 1.5.2      2025-09-08 [1] CRAN (R 4.3.0)
##  SummarizedExperiment * 1.32.0     2023-10-24 [1] Bioconductor
##  svglite                2.1.3      2023-12-08 [1] CRAN (R 4.3.0)
##  systemfonts            1.2.1      2025-01-20 [1] CRAN (R 4.3.0)
##  tibble               * 3.3.0      2025-06-08 [1] CRAN (R 4.3.0)
##  tidygraph              1.3.1      2024-01-30 [1] CRAN (R 4.3.0)
##  tidyr                * 1.3.1      2024-01-24 [1] CRAN (R 4.3.0)
##  tidyselect             1.2.1      2024-03-11 [1] CRAN (R 4.3.0)
##  tidytree               0.4.6      2023-12-12 [1] CRAN (R 4.3.0)
##  tidyverse            * 2.0.0      2023-02-22 [1] CRAN (R 4.3.0)
##  timechange             0.3.0      2024-01-18 [1] CRAN (R 4.3.0)
##  treeio                 1.26.0     2023-10-24 [1] Bioconductor
##  tweenr                 2.0.3      2024-02-26 [1] CRAN (R 4.3.0)
##  tzdb                   0.4.0      2023-05-12 [1] CRAN (R 4.3.0)
##  variancePartition    * 1.32.5     2024-02-16 [1] Bioconductor 3.18 (R 4.3.0)
##  vctrs                  0.6.5      2023-12-01 [1] CRAN (R 4.3.0)
##  viridis                0.6.5      2024-01-29 [1] CRAN (R 4.3.0)
##  viridisLite            0.4.2      2023-05-02 [1] CRAN (R 4.3.0)
##  vroom                  1.6.5      2023-12-05 [1] CRAN (R 4.3.0)
##  withr                  3.0.2      2024-10-28 [1] CRAN (R 4.3.0)
##  xfun                   0.53       2025-08-19 [1] CRAN (R 4.3.0)
##  xml2                   1.4.0      2025-08-20 [1] CRAN (R 4.3.0)
##  XVector                0.42.0     2023-10-24 [1] Bioconductor
##  yaml                   2.3.10     2024-07-26 [1] CRAN (R 4.3.0)
##  yulab.utils            0.1.9      2025-01-07 [1] CRAN (R 4.3.0)
##  zlibbioc               1.48.2     2024-03-13 [1] Bioconductor 3.18 (R 4.3.0)
## 
##  [1] /opt/R/4.3.0/lib/R/library
## 
## ──────────────────────────────────────────────────────────────────────────────

PD Differential Gene Expression - DEG with Limma/voom Hardy project dataset

Corrected identification script

Guillermo Rocamora Pérez

30 September 2025

1 JB changes

2 Background

2.1 Load previous data

Changelog

v1.6

v1.5

v1.4

v1.3

v1.2

v1.1

v1.0

3 DGE with limma/voom and duplicateCorrelation

3.1 Preprocessing

3.2 Voom with duplicateCorrelation

3.3 Fit model

3.4 Make constrasts

4 Results: Collapsed brain region

4.1 Volcano plots

PD vs Control

4.2 Top significant genes

PD vs Control

4.3 Mean-difference plots

PD vs Control

5 Results: By brain region

5.1 Volcano plots

ACG

IPL

MFG

MTG

5.2 Top Significant genes

ACG

IPL

MFG

MTG

5.3 Mean-difference plots

ACG

IPL

MFG

MTG

6 Summary of results

6.1 Number of differentially expressed genes

PD vs Control

6.2 Common DEGs between brain regions

7 Output results to csv

8 Session Info

PD Differential Gene Expression - DEG with Limma/voom
Hardy project dataset

3 DGE with `limma/voom` and `duplicateCorrelation`

3.2 Voom with `duplicateCorrelation`